mixlink (version 0.1.5)

Metropolis-Hastings Sampler: Metropolis-Hastings Sampler

Description

Random-Walk Metropolis Hastings sampler for Binomial and Poisson Mixture Link models.

Usage

rwmetrop.mixlink.binomial(y, X, m, R, burn = 0, thin = 1,
  invlink = plogis, Beta.init = NULL, Pi.init, kappa.init = NULL,
  hyper = NULL, report.period = R + 1, use.laplace.approx = TRUE,
  proposal = NULL, param.grp = NULL, fixed.kappa = FALSE)
rwmetrop.mixlink.poisson(y, X, offset = rep(0, length(y)), R, burn = 0,
  thin = 1, invlink = exp, Beta.init = NULL, Pi.init, kappa.init = NULL,
  hyper = NULL, report.period = R + 1, use.laplace.approx = TRUE,
  proposal = NULL, param.grp = NULL, fixed.kappa = FALSE)

Arguments

y

Observations.

X

Design matrix for regression.

m

Number of success/failure trials.

R

Number of MCMC draws to take.

burn

Number of initial MCMC draws to discard.

thin

After burn-in period, save one of every thin draws.

invlink

The inverse link function for the mean. Default is plogis for Binomial and exp for Poisson.

Beta.init

Starting value for \(\bm{\beta}\).

Pi.init

Starting value for \(\bm{\pi}\).

kappa.init

Starting value for \(\kappa\).

hyper

A list with hyperparameters corresponding to the prior from the Details section. var.Beta is \(V_{\beta}\) with default 1000. alpha.Pi is \(\bm{\alpha}_\pi\) with default rep(1,J). a.kappa and b.kappa correspond to \((a_\kappa, b_\kappa)\) which have default values a.kappa = 1 and b.kappa = 1/10.

report.period

Report progress every report.period draws.

use.laplace.approx

Maximize a Laplace approximation to the posterior, to find a starting value for MCMC.

proposal

A list with two elements. var is the covariance matrix for a \(d + (J-1) + 1\) dimensional multivariate normal proposal, scale is a scalar multiplied with var. Defaults are var = diag(d+(J-1)+1) and scale = 0.02.

param.grp

A vector of integers of length \(d + (J-1) + 1\), where d = ncol(X), which indicates the grouping of parameters in MCMC. Parameters with common integers are sampled together. The first \(d\) correspond to \(\bm{\beta}\), the next \(J-1\) correspond to \(\bm{\pi}\), and the last one corresponds to \(\kappa\). At the default value (param.grp = NULL), \(\bm{\beta}\), \(\bm{\pi}\), and \(\kappa\) are sampled one at a time, each in their entirety.

fixed.kappa

Keep \(\kappa\) fixed at kappa.init (Default FALSE).

offset

Constant offset term to add to \(x^T \beta\).

Value

A list with the MCMC results:

par.hist

\(R \times [d + (J-1) + 1]\) matrix of saved MCMC draws before transformations are applied. Most users will not need this.

Beta.hist

\(R \times d\) matrix of saved \(\bm{\beta}\) draws

Pi.hist

\(R \times J\) matrix of saved \(\bm{\pi}\) draws

kappa.hist

\(R \times 1\) vector of \(\kappa\) draws

accept

Percentages that MCMC proposals were accepted. Corresponds to param.grp

elapsed.sec

Elapsed time for sampling, in seconds.

laplace.out

Output of Laplace approximation.

R.keep

Number of draws kept, after thinning and burn-in.

X.names

Names of columns of design matrix.

Can be accessed with the functions print, summary, and DIC.

Details

Priors for Bayesian Mixture Link model are

• \( \bm{\beta} \sim \textrm{N}(\bm{0}, V_{\beta} \bm{I}) \),

• \( \bm{\pi} \sim \textrm{Dirichlet}_J(\bm{\alpha}_\pi) \),

• \( \kappa \sim \textrm{Gamma}(a_\kappa, b_\kappa) \), parameterized with \(\textrm{E}(\kappa) = a_\kappa / b_\kappa\).

References

Andrew M. Raim, Nagaraj K. Neerchal, and Jorge G. Morel. An Extension of Generalized Linear Models to Finite Mixture Outcomes. arXiv preprint: 1612.03302

Examples

# NOT RUN {
library(Matrix)

# ----- Generate data -----
n <- 200
x <- runif(n, 1, 3)
X <- model.matrix(~ x)
Beta.true <- c(0, 1)
mean.true <- exp(X %*% Beta.true)
kappa.true <- 0.95
Pi.true <- c(1,3)/4
d <- ncol(X)
J <- length(Pi.true)
y <- r.mixlink.pois(n, mean.true, Pi.true, kappa.true)

# ----- Run Metropolis-within-Gibbs sampler -----
hyper <- list(VBeta = diag(1000, d), alpha.Pi = rep(1, J),
	kappa.a = 1, kappa.b = 1/2)
proposal <- list(
var = bdiag(solve(t(X) %*% X), diag(J-1), 1),
scale = 0.5)
metrop.out <- rwmetrop.mixlink.poisson(y, X, R = 20000, burn = 1000,
	thin = 10, Pi.init = c(1,9)/10, hyper = hyper,
	report.period = 1000, use.laplace.approx = TRUE, proposal = proposal)

print(metrop.out)
DIC(metrop.out)
# }
# NOT RUN {
# }